#######################################################################################
### R file to replicate the creation of the figures found in "Long-Term Effects in Models with Temporal Dependence", Political Analysis, Laron K. Williams
###
### NOTE: email williamslaro@missouri.edu if you have any questions.
#######################################################################################

library(foreign)
library(ggplot2)
library(lattice)
library(fields)
library(separationplot)

### Set up working directory containing all of the files.
#setwd("")

#######################################################################################
### Figure 1: Variations of long-term effects based on the counterfactual of interest for negative duration dependence
#######################################################################################
lte <- read.dta("LTE Example--Negative.dta", convert.underscore = TRUE)

### 1-time shift
n1 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 13), aes(x=xaxis, y=pr.s1.0.0)) + xlab("") + ylab("Pr(Y=1)")
n1 <- n1 + geom_point() + geom_linerange(aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo), linetype="dashed")
n1 <- n1 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s1.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s1.1.0.hi, ymin = pr.s1.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n1 <- n1 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo, linetype="Y0"))
n1 <- n1 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n1 <- n1 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s1.1.0), vjust=-0.25, hjust=-.25, size=3.5)
n1 <- n1 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n1 <- n1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n1 <- n1 + scale_y_continuous(limits=c(-0.01, 0.4), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n1

### Temporary (4-period) shift
n2 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s2.0.0)) + xlab("") + ylab("Pr(Y=1)")
n2 <- n2 + geom_point() + geom_linerange(aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo), linetype="dashed")
n2 <- n2 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s2.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s2.1.0.hi, ymin = pr.s2.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n2 <- n2 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo, linetype="Y0"))
n2 <- n2 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n2 <- n2 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s2.1.0), vjust=-0.25, hjust=-.25, size=3.5)
n2 <- n2 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n2 <- n2 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n2 <- n2 + scale_y_continuous(limits=c(-0.01, 0.4), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n2

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n3 <- n3 + scale_y_continuous(limits=c(-0.01, 0.4), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n3

### Compounded effects
n4 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n4 <- n4 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n4 <- n4 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n4 <- n4 + geom_point(data=subset(lte, xaxis > 5), aes(y=pr.s3.1.1, x=xaxis-0.3)) + geom_linerange(data=subset(lte, xaxis > 5), aes(ymax = pr.s3.1.1.hi, ymin = pr.s3.1.1.lo, x=xaxis-0.3, linetype="Y11"))
n4 <- n4 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=1, hjust=-.5, size=3.5)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 5), aes(label = time1.1, y=pr.s3.1.1), vjust=-0.75, hjust=3.5, size=3.5)
n4 <- n4 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n4 <- n4 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1","Y11"), labels=c(expression(Y[t]==0), expression(Y[t]==1), expression(Y[t+5]==1)), values=c("dashed","solid","dotted")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n4 <- n4 + scale_y_continuous(limits=c(-0.01, 0.4), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n4


#######################################################################################
### Figure 2a: Long-term effects with non-proportional hazards ignored
#######################################################################################
lte <- read.dta("LTE Example--Negative with Positive NPH.dta", convert.underscore = TRUE)

n3_n <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 17), aes(x=xaxis, y=pr.0.0.n.mn)) + xlab("") + ylab("Pr(Y=1)")
n3_n <- n3_n + geom_point() + geom_linerange(aes(ymax = pr.0.0.n.hi, ymin = pr.0.0.n.lo), linetype="dashed")
n3_n <- n3_n + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.0.1.n.mn, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.0.1.n.hi, ymin = pr.0.1.n.lo, x=xaxis-0.1, linetype="Y1"))
n3_n <- n3_n + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.0.0.n.hi, ymin = pr.0.0.n.lo, linetype="Y0"))
n3_n <- n3_n + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n3_n <- n3_n + geom_text(data=subset(lte, xaxis > 0), aes(label = time0.1, y=pr.0.1.n.mn), vjust=-0.25, hjust=-.25, size=3.5)
n3_n <- n3_n + scale_x_continuous(limits=c(-2,16), breaks=c(-1, 0, 4, 8, 12, 16), labels=c("Baseline", "t", "t+4", "t+8", "t+12", "t+16"))
n3_n <- n3_n + scale_y_continuous(limits=c(0,.45), breaks=c(0, .1, .2, .3, .4))
n3_n <- n3_n + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n3_n


#######################################################################################
### Figure 2b: Long-term effects with non-proportional hazards modeled
#######################################################################################
n3_y <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 17), aes(x=xaxis, y=pr.0.0.y.mn)) + xlab("") + ylab("Pr(Y=1)")
n3_y <- n3_y + geom_point() + geom_linerange(aes(ymax = pr.0.0.y.hi, ymin = pr.0.0.y.lo), linetype="dashed")
n3_y <- n3_y + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.0.1.y.mn, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.0.1.y.hi, ymin = pr.0.1.y.lo, x=xaxis-0.1, linetype="Y1"))
n3_y <- n3_y + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.0.0.y.hi, ymin = pr.0.0.y.lo, linetype="Y0"))
n3_y <- n3_y + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n3_y <- n3_y + geom_text(data=subset(lte, xaxis > 0), aes(label = time0.1, y=pr.0.1.y.mn), vjust=-0.25, hjust=-.25, size=3.5)
n3_y <- n3_y + scale_x_continuous(limits=c(-2,16), breaks=c(-1, 0, 4, 8, 12, 16), labels=c("Baseline", "t", "t+4", "t+8", "t+12", "t+16"))
n3_y <- n3_y + scale_y_continuous(limits=c(0,.45), breaks=c(0, .1, .2, .3, .4))
n3_y <- n3_y + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75), legend.text=element_text(size=14), axis.text=element_text(size=12))
n3_y


#######################################################################################
### Figure 2c: Difference in Long-Term Effects: (LTE | NPH Modeled) - (LTE | NPH Ignored)
#######################################################################################
n3 <- ggplot(subset(lte, xaxis >= 1 & xaxis <= 17), aes(x=xaxis, y=diff.lte.mn)) + xlab("") + ylab("(LTE | NPH Modeled) - (LTE | NPH Ignored)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = diff.lte.hi, ymin = diff.lte.lo), linetype="solid")
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3.5)
n3 <- n3 + scale_x_continuous(limits=c(0,16), breaks=c(1, 4, 8, 12, 16), labels=c("t+1", "t+4", "t+8", "t+12", "t+16"))
n3 <- n3 + geom_hline(y=0, linetype="dashed") + theme(axis.text=element_text(size=12))
n3


#######################################################################################
### Figure 3: Average hazard rates compared to true hazard rates
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_1.dta", convert.underscore = TRUE)
mn_nonph <- subset(mn, scenario != "Negative with Positive NPH" & scenario != "Positive with Negative NPH")

p <- ggplot(subset(mn_nonph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85), axis.text=element_text(size=12), strip.text.x = element_text(size = 16))
p <- p + scale_x_continuous(limits=c(0,25), breaks=c(0, 5, 10, 15, 20))
p



#######################################################################################
### Figure 4: Probabilities of long-term effects when ideological distance is at its minimum and its maximum values
#######################################################################################
ispp = read.dta("ISPP.dta", convert.underscore=TRUE)

### Minimum value of ideological distance
t <- ggplot(ispp, aes(t, p.min)) + geom_point(alpha = 0.2) + ylab("Pr(Initiation)") + xlab("Time Since Previous Dispute")
t <- t + theme(axis.text=element_text(size=12))
t <- t + ylim(0, 0.3) #+ scale_x_continuous(breaks = seq(0, 3000, 1000))
t

### Maximum value of ideological distance
t2 <- ggplot(ispp, aes(t, p.max)) + geom_point(alpha = 0.2) + ylab("Pr(Initiation)") + xlab("Time Since Previous Dispute")
t2 <- t2 + theme(axis.text=element_text(size=12))
t2 <- t2 + ylim(0, 0.3) #+ scale_x_continuous(breaks = seq(0, 3000, 1000))
t2


#######################################################################################
### Figure 5: Average long-term effects across time since previous dispute
#######################################################################################
alte = read.dta("Clare--LTE--Average.dta", convert.underscore=TRUE)

alte$xaxis <- alte$xaxis - 1
alte$lab <- alte$xaxis - 1
alte$time0 <- alte$time0 - 1
alte$time1 <- alte$time1 - 1

# Levels
shift <- 0.2
a1 <- ggplot(subset(alte, xaxis >= 0 & xaxis <= 20), aes(x=xaxis, y=p0)) + xlab("") + ylab("Pr(Y=1)")
a1 <- a1 + geom_point() + geom_linerange(aes(ymax = p0.hi, ymin = p0.lo), linetype="dashed")
a1 <- a1 + geom_point(data=subset(alte, xaxis > 0), aes(y=p1, x=xaxis-shift)) + geom_linerange(data=subset(alte, xaxis > 0), aes(ymax = p1.hi, ymin = p1.lo, x=xaxis-shift, linetype="Y1"))
a1 <- a1 + geom_point(data=subset(alte, xaxis < 1)) + geom_linerange(data=subset(alte, xaxis < 1), aes(ymax = p0.hi, ymin = p0.lo, linetype="Y0"))
a1 <- a1 + geom_text(aes(label = time0), vjust=-0.25, hjust=-0.25, size=3.25)
a1 <- a1 + geom_text(data=subset(alte, xaxis > 0), aes(label = time1, y= p1), vjust=-0.25, hjust=2, size=3)
a1 <- a1 + scale_x_continuous(limits=c(0,20), breaks=c(0, 5, 10, 15, 20), labels=c("t", "t+5", "t+10", "t+15", "t+20"))
a1 <- a1 + scale_y_continuous(limits=c(-0.01,.31), breaks=c(0, .1, .2, .3))
a1 <- a1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75))
a1 <- a1 + theme(axis.text=element_text(size=12))
a1

# Difference
cdci <- ggplot(data=subset(alte, xaxis >= 1 & xaxis <= 20), aes(x=xaxis)) + geom_point(aes(y = lte), size=1.5) + geom_linerange(aes(ymax = lte.hi, ymin = lte.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed")
cdci <- cdci + geom_text(aes(y=lte, label = lab), vjust=0, hjust=-0.5, size=3.25)
cdci <- cdci + xlab("") + ylab("Long-Term Effect")
cdci <- cdci + scale_x_continuous(limits=c(0,20), breaks=c(0, 5, 10, 15, 20), labels=c("t", "t+5", "t+10", "t+15", "t+20"))
cdci <- cdci + scale_y_continuous(limits=c(-0.01,.31), breaks=c(0, .1, .2, .3))
cdci <- cdci + theme(axis.text=element_text(size=12))
cdci




#######################################################################################
### Figure 6: Long-term effects across time since previous dispute for four simulation scenarios with varying baseline probabilities
#######################################################################################
slte = read.dta("Clare--LTE--Scenarios.dta", convert.underscore=TRUE)

slte$xaxis <- slte$xaxis - 1
slte$time1 <- slte$time1 - 1

slte1 <- subset(slte, scenario == "no")
slte2 <- subset(slte, scenario == "low")
slte3 <- subset(slte, scenario == "mod")
slte4 <- subset(slte, scenario == "hi")

hj <- -0.5
vj <- -0.25
sz <- 3.5

### Panel a: Pr = 0.003
p1 <- ggplot(data=subset(slte1, xaxis > 0), aes(x = xaxis)) + geom_point(aes(y = lte.s1), size = 1) + xlab("") + ylab("Long-Term Effect") 
p1 <- p1 + geom_linerange(aes(ymin=lte.s1.lo, ymax=lte.s1.hi), type="solid")
p1 <- p1 + geom_text(aes(y=lte.s1, label = time1), vjust=vj, hjust=hj, size=sz)
p1 <- p1 + geom_hline(yintercept = 0, linetype = "dashed") + theme(axis.text=element_text(size=12))
p1 <- p1 + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
p1 <- p1 + ylim(-0.005, 0.2)
p1

### Panel b: Pr = 0.05
p2 <- ggplot(data=subset(slte2, xaxis > 0), aes(x = xaxis)) + geom_point(aes(y = lte.s1), size = 1) + xlab("") + ylab("Long-Term Effect") 
p2 <- p2 + geom_linerange(aes(ymin=lte.s1.lo, ymax=lte.s1.hi), type="solid")
p2 <- p2 + geom_text(aes(y=lte.s1, label = time1), vjust=vj, hjust=hj, size=sz)
p2 <- p2 + geom_hline(yintercept = 0, linetype = "dashed") + theme(axis.text=element_text(size=12))
p2 <- p2 + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
p2 <- p2 + ylim(-0.005, 0.2)
p2

### Panel c: Pr = 0.10
p3 <- ggplot(data=subset(slte3, xaxis > 0), aes(x = xaxis)) + geom_point(aes(y = lte.s1), size = 1) + xlab("") + ylab("Long-Term Effect") 
p3 <- p3 + geom_linerange(aes(ymin=lte.s1.lo, ymax=lte.s1.hi), type="solid")
p3 <- p3 + geom_text(aes(y=lte.s1, label = time1), vjust=vj, hjust=hj, size=sz)
p3 <- p3 + geom_hline(yintercept = 0, linetype = "dashed") + theme(axis.text=element_text(size=12))
p3 <- p3 + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
p3 <- p3 + ylim(-0.005, 0.2)
p3

### Panel d: Pr = 0.50
p4 <- ggplot(data=subset(slte4, xaxis > 0), aes(x = xaxis)) + geom_point(aes(y = lte.s1), size = 1) + xlab("") + ylab("Long-Term Effect") 
p4 <- p4 + geom_linerange(aes(ymin=lte.s1.lo, ymax=lte.s1.hi), type="solid")
p4 <- p4 + geom_text(aes(y=lte.s1, label = time1), vjust=vj, hjust=hj, size=sz)
p4 <- p4 + geom_hline(yintercept = 0, linetype = "dashed") + theme(axis.text=element_text(size=12))
p4 <- p4 + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
p4 <- p4 + ylim(-0.005, 0.2)
p4





#######################################################################################
#######################################################################################
########################### APPENDIX FIGURES ##########################################
#######################################################################################
#######################################################################################

#######################################################################################
### Figure S.1: Four functional forms of duration dependence
#######################################################################################
ff = read.dta("Functional Forms.dta", convert.underscore = TRUE)

p <- ggplot(ff, aes(x=h.t, y=h, group=scen)) + geom_line(aes(linetype=scen))
p <- p + xlab("t") + ylab("Pr(y=1)") + annotate("text", label = "Positive", x = 2.5, y = 0.07, size=4) + annotate("text", label = "Negative", x = 22, y = 0.085, size=4) + annotate("text", label = "Non-Monotonic 1", x = 3.75, y = 0.32, size=4) + annotate("text", label = "Non-Monotonic 2", x = 15.5, y = 0.2, size=4)
p <- p + scale_linetype_manual(values=c(2, 1, 3, 4))
p <- p + theme(legend.position="none")
p


#######################################################################################
### Figure S.2: Long-term effects for four scenarios for negative duration dependence
#######################################################################################
ndd = read.dta("LTE Example--Scenarios--Negative.dta", convert.underscore = TRUE)

sc1 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y=lte.s1), size=1.5) + geom_linerange(aes(ymax = lte.s1.hi, ymin = lte.s1.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc1 <- sc1 + geom_text(aes(y=lte.s1, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc1 <- sc1 + ggtitle(expression(paste(t["XC"], "=1"))) + xlim(0, 20) + ylim(-0.05, 0.25)
sc1

sc2 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y = lte.s2), size=1.5) + geom_linerange(aes(ymax = lte.s2.hi, ymin = lte.s2.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc2 <- sc2 + geom_text(aes(y=lte.s2, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc2 <- sc2 + ggtitle(expression(paste(t["XC"], "=5"))) + xlim(0, 20) + ylim(-0.05, 0.25)
sc2

sc3 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 15), aes(x=xaxis)) + geom_point(aes(y = lte.s3), size=1.5) + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc3 <- sc3 + geom_text(aes(y=lte.s3, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc3 <- sc3 + ggtitle(expression(paste(t["XC"], "=10"))) + xlim(0, 20) + ylim(-0.05, 0.25) 
sc3

sc4 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 10), aes(x=xaxis)) + geom_point(aes(y = lte.s4), size=1.5) + geom_linerange(aes(ymax = lte.s4.hi, ymin = lte.s4.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc4 <- sc4 + geom_text(aes(y=lte.s4, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc4 <- sc4 + ggtitle(expression(paste(t["XC"], "=20"))) + xlim(0, 20) + ylim(-0.05, 0.25) 
sc4

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(sc1, vp = vplayout(1, 1))
print(sc2, vp = vplayout(1, 2))
print(sc3, vp = vplayout(2, 1))
print(sc4, vp = vplayout(2, 2))

#######################################################################################
### Figure S.3: Long-term effects for four scenarios for positive duration dependence
#######################################################################################
ndd = read.dta("LTE Example--Scenarios--Positive.dta", convert.underscore = TRUE)

sc1 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y=lte.s1), size=1.5) + geom_linerange(aes(ymax = lte.s1.hi, ymin = lte.s1.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc1 <- sc1 + geom_text(aes(y=lte.s1, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc1 <- sc1 + ggtitle(expression(paste(t["XC"], "=1"))) + xlim(0, 20) + ylim(-1, 0.05)
sc1

sc2 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y = lte.s2), size=1.5) + geom_linerange(aes(ymax = lte.s2.hi, ymin = lte.s2.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc2 <- sc2 + geom_text(aes(y=lte.s2, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc2 <- sc2 + ggtitle(expression(paste(t["XC"], "=5"))) + xlim(0, 20) + ylim(-1, 0.05)
sc2

sc3 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 15), aes(x=xaxis)) + geom_point(aes(y = lte.s3), size=1.5) + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc3 <- sc3 + geom_text(aes(y=lte.s3, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc3 <- sc3 + ggtitle(expression(paste(t["XC"], "=10"))) + xlim(0, 20) + ylim(-1, 0.05)
sc3

sc4 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 10), aes(x=xaxis)) + geom_point(aes(y = lte.s4), size=1.5) + geom_linerange(aes(ymax = lte.s4.hi, ymin = lte.s4.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc4 <- sc4 + geom_text(aes(y=lte.s4, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc4 <- sc4 + ggtitle(expression(paste(t["XC"], "=20"))) + xlim(0, 20) + ylim(-1, 0.05)
sc4

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(sc1, vp = vplayout(1, 1))
print(sc2, vp = vplayout(1, 2))
print(sc3, vp = vplayout(2, 1))
print(sc4, vp = vplayout(2, 2))



#######################################################################################
### Figure S.4: Long-term effects for four scenarios for non-monotonic (1) duration dependence
#######################################################################################
ndd = read.dta("LTE Example--Scenarios--NM1.dta", convert.underscore = TRUE)

sc1 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y=lte.s1), size=1.5) + geom_linerange(aes(ymax = lte.s1.hi, ymin = lte.s1.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc1 <- sc1 + geom_text(aes(y=lte.s1, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc1 <- sc1 + ggtitle(expression(paste(t["XC"], "=1"))) + xlim(0, 20) + ylim(-0.5, 0.5)
sc1

sc2 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y = lte.s2), size=1.5) + geom_linerange(aes(ymax = lte.s2.hi, ymin = lte.s2.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc2 <- sc2 + geom_text(aes(y=lte.s2, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc2 <- sc2 + ggtitle(expression(paste(t["XC"], "=5"))) + xlim(0, 20) + ylim(-0.5, 0.5)
sc2

sc3 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 15), aes(x=xaxis)) + geom_point(aes(y = lte.s3), size=1.5) + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc3 <- sc3 + geom_text(aes(y=lte.s3, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc3 <- sc3 + ggtitle(expression(paste(t["XC"], "=10"))) + xlim(0, 20) + ylim(-0.5, 0.5)
sc3

sc4 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 10), aes(x=xaxis)) + geom_point(aes(y = lte.s4), size=1.5) + geom_linerange(aes(ymax = lte.s4.hi, ymin = lte.s4.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc4 <- sc4 + geom_text(aes(y=lte.s4, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc4 <- sc4 + ggtitle(expression(paste(t["XC"], "=20"))) + xlim(0, 20) + ylim(-0.5, 0.5)
sc4

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(sc1, vp = vplayout(1, 1))
print(sc2, vp = vplayout(1, 2))
print(sc3, vp = vplayout(2, 1))
print(sc4, vp = vplayout(2, 2))


#######################################################################################
### Figure S.5: Long-term effects for four scenarios for non-monotonic (2) duration dependence
#######################################################################################
ndd = read.dta("LTE Example--Scenarios--NM2.dta", convert.underscore = TRUE)

sc1 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y=lte.s1), size=1.5) + geom_linerange(aes(ymax = lte.s1.hi, ymin = lte.s1.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc1 <- sc1 + geom_text(aes(y=lte.s1, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc1 <- sc1 + ggtitle(expression(paste(t["XC"], "=1"))) + xlim(0, 20) + ylim(-0.2, 0.3)
sc1

sc2 <- ggplot(data=subset(ndd, xaxis > 0), aes(x=xaxis)) + geom_point(aes(y = lte.s2), size=1.5) + geom_linerange(aes(ymax = lte.s2.hi, ymin = lte.s2.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc2 <- sc2 + geom_text(aes(y=lte.s2, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc2 <- sc2 + ggtitle(expression(paste(t["XC"], "=5"))) + xlim(0, 20) + ylim(-0.2, 0.3)
sc2

sc3 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 15), aes(x=xaxis)) + geom_point(aes(y = lte.s3), size=1.5) + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo), linetype="solid") + xlab("") + ylab("Long-Term Effect") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc3 <- sc3 + geom_text(aes(y=lte.s3, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc3 <- sc3 + ggtitle(expression(paste(t["XC"], "=10"))) + xlim(0, 20) + ylim(-0.2, 0.3)
sc3

sc4 <- ggplot(data=subset(ndd, xaxis > 0 & xaxis <= 10), aes(x=xaxis)) + geom_point(aes(y = lte.s4), size=1.5) + geom_linerange(aes(ymax = lte.s4.hi, ymin = lte.s4.lo), linetype="solid") + xlab("") + ylab("") + geom_hline(yintercept = 0, linetype = "dashed", colour = "red")
sc4 <- sc4 + geom_text(aes(y=lte.s4, label = time1), vjust=-0.75, hjust=-0.75, size=2)
sc4 <- sc4 + ggtitle(expression(paste(t["XC"], "=20"))) + xlim(0, 20) + ylim(-0.2, 0.3)
sc4


grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(sc1, vp = vplayout(1, 1))
print(sc2, vp = vplayout(1, 2))
print(sc3, vp = vplayout(2, 1))
print(sc4, vp = vplayout(2, 2))

#######################################################################################
### Figure S.6: Two methods of depicting long-term effects: negative duration dependence
#######################################################################################
lte <- read.dta("LTE Example--Negative.dta", convert.underscore = TRUE)

n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75))
n3 <- n3 + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n3

### LTE: permanent effects
lten <- ggplot(subset(lte, xaxis >=1 & xaxis <= 12), aes(x=xaxis, y=lte.s3)) + xlab("") + ylab("Pr(Y=1)")
lten <- lten + geom_point() + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo))
lten <- lten + geom_point(data=subset(lte, xaxis == 0), aes(y=pr.s3.1.0)) + geom_linerange(data=subset(lte, xaxis == 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo))
lten <- lten + geom_hline(y=0, linetype=2)
lten <- lten + geom_text(aes(label = time1.0, y=lte.s3), vjust=-0.25, hjust=-1, size=3)
lten <- lten + geom_text(data=subset(lte, xaxis == 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-1, size=3)
lten <- lten + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
lten <- lten + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
lten



#######################################################################################
### Figure S.7: Two methods of depicting long-term effects: positive duration dependence
#######################################################################################
lte <- read.dta("LTE Example--Positive.dta", convert.underscore = TRUE)

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.25, 0.75))
n3

### LTE: permanent effects
ltep <- ggplot(subset(lte, xaxis >=1 & xaxis <= 12), aes(x=xaxis, y=lte.s3)) + xlab("") + ylab("Pr(Y=1)")
ltep <- ltep + geom_point() + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo))
ltep <- ltep + geom_point(data=subset(lte, xaxis == 0), aes(y=pr.s3.1.0)) + geom_linerange(data=subset(lte, xaxis == 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo))
ltep <- ltep + geom_hline(y=0, linetype=2)
ltep <- ltep + geom_text(aes(label = time1.0, y=lte.s3), vjust=-0.25, hjust=-0.5, size=3)
ltep <- ltep + geom_text(data=subset(lte, xaxis == 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-1, size=3)
ltep <- ltep + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
ltep

#######################################################################################
### Figure S.8: Two methods of depicting long-term effects: non-monotonic duration dependence (parabolic)
#######################################################################################
lte <- read.dta("LTE Example--NM1.dta", convert.underscore = TRUE)

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.75, 0.75))
n3

### LTE: permanent effects
ltep <- ggplot(subset(lte, xaxis >=1 & xaxis <= 12), aes(x=xaxis, y=lte.s3)) + xlab("") + ylab("Pr(Y=1)")
ltep <- ltep + geom_point() + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo))
ltep <- ltep + geom_point(data=subset(lte, xaxis == 0), aes(y=pr.s3.1.0)) + geom_linerange(data=subset(lte, xaxis == 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo))
ltep <- ltep + geom_hline(y=0, linetype=2)
ltep <- ltep + geom_text(aes(label = time1.0, y=lte.s3), vjust=-0.25, hjust=-0.5, size=3)
ltep <- ltep + geom_text(data=subset(lte, xaxis == 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-1, size=3)
ltep <- ltep + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
ltep


#######################################################################################
### Figure S.9: Two methods of depicting long-term effects: non-monotonic duration dependence (log-logistic)
#######################################################################################
lte <- read.dta("LTE Example--NM2.dta", convert.underscore = TRUE)

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.9, 0.9))
n3

### LTE: permanent effects
ltep <- ggplot(subset(lte, xaxis >=1 & xaxis <= 12), aes(x=xaxis, y=lte.s3)) + xlab("") + ylab("Pr(Y=1)")
ltep <- ltep + geom_point() + geom_linerange(aes(ymax = lte.s3.hi, ymin = lte.s3.lo))
ltep <- ltep + geom_point(data=subset(lte, xaxis == 0), aes(y=pr.s3.1.0)) + geom_linerange(data=subset(lte, xaxis == 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo))
ltep <- ltep + geom_hline(y=0, linetype=2)
ltep <- ltep + geom_text(aes(label = time1.0, y=lte.s3), vjust=-0.25, hjust=-0.5, size=3)
ltep <- ltep + geom_text(data=subset(lte, xaxis == 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-1, size=3)
ltep <- ltep + scale_x_continuous(limits=c(0,12), breaks=c(0, 4, 8, 12), labels=c("t", "t+4", "t+8", "t+12"))
ltep


#######################################################################################
### Figure S.10: Variations of long-term effects based on the counterfactual of interest for negative duration dependence
#######################################################################################
lte <- read.dta("LTE Example--Negative.dta", convert.underscore = TRUE)

### 1-time shift in effects
n1 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 13), aes(x=xaxis, y=pr.s1.0.0)) + xlab("") + ylab("Pr(Y=1)")
n1 <- n1 + geom_point() + geom_linerange(aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo), linetype="dashed")
n1 <- n1 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s1.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s1.1.0.hi, ymin = pr.s1.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n1 <- n1 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo, linetype="Y0"))
n1 <- n1 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n1 <- n1 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s1.1.0), vjust=-0.25, hjust=-.25, size=3)
n1 <- n1 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n1 <- n1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75))
n1 <- n1 + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n1

### Temporary (4-period) shift in effects
n2 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s2.0.0)) + xlab("") + ylab("Pr(Y=1)")
n2 <- n2 + geom_point() + geom_linerange(aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo), linetype="dashed")
n2 <- n2 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s2.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s2.1.0.hi, ymin = pr.s2.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n2 <- n2 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo, linetype="Y0"))
n2 <- n2 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n2 <- n2 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s2.1.0), vjust=-0.25, hjust=-.25, size=3)
n2 <- n2 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n2 <- n2 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75))
n2 <- n2 + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n2

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.85, 0.75))
n3 <- n3 + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n3

### Compounded effects
n4 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n4 <- n4 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n4 <- n4 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n4 <- n4 + geom_point(data=subset(lte, xaxis > 5), aes(y=pr.s3.1.1, x=xaxis-0.3)) + geom_linerange(data=subset(lte, xaxis > 5), aes(ymax = pr.s3.1.1.hi, ymin = pr.s3.1.1.lo, x=xaxis-0.3, linetype="Y11"))
n4 <- n4 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=1, hjust=-.5, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 5), aes(label = time1.1, y=pr.s3.1.1), vjust=-0.75, hjust=3.5, size=3)
n4 <- n4 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n4 <- n4 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1","Y11"), labels=c(expression(Y[t]==0), expression(Y[t]==1), expression(Y[t+5]==1)), values=c("dashed","solid","dotted")) + theme(legend.position=c(0.85, 0.75))
n4 <- n4 + scale_y_continuous(limits=c(-0.01, 0.475), breaks=c(0, 0.1, 0.2, 0.3, 0.4))
n4



#######################################################################################
### Figure S.11: Variations of long-term effects based on the counterfactual of interest for positive duration dependence
#######################################################################################
lte <- read.dta("LTE Example--Positive.dta", convert.underscore = TRUE)

### 1-time shift in effects
n1 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s1.0.0)) + xlab("") + ylab("Pr(Y=1)")
n1 <- n1 + geom_point() + geom_linerange(aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo), linetype="dashed")
n1 <- n1 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s1.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s1.1.0.hi, ymin = pr.s1.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n1 <- n1 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo, linetype="Y0"))
n1 <- n1 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n1 <- n1 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s1.1.0), vjust=-0.25, hjust=-.25, size=3)
n1 <- n1 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n1 <- n1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.25, 0.75))
n1

### Temporary (4-period) shift in effects
n2 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s2.0.0)) + xlab("") + ylab("Pr(Y=1)")
n2 <- n2 + geom_point() + geom_linerange(aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo), linetype="dashed")
n2 <- n2 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s2.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s2.1.0.hi, ymin = pr.s2.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n2 <- n2 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo, linetype="Y0"))
n2 <- n2 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n2 <- n2 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s2.1.0), vjust=-0.25, hjust=-.25, size=3)
n2 <- n2 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n2 <- n2 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.25, 0.75))
n2

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.25, 0.75))
n3

### Compounded effects
n4 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n4 <- n4 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n4 <- n4 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n4 <- n4 + geom_point(data=subset(lte, xaxis > 5), aes(y=pr.s3.1.1, x=xaxis-0.3)) + geom_linerange(data=subset(lte, xaxis > 5), aes(ymax = pr.s3.1.1.hi, ymin = pr.s3.1.1.lo, x=xaxis-0.3, linetype="Y11"))
n4 <- n4 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.75, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 5), aes(label = time1.1, y=pr.s3.1.1), vjust=1, hjust=3.5, size=3)
n4 <- n4 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n4 <- n4 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1","Y11"), labels=c(expression(Y[t]==0), expression(Y[t]==1), expression(Y[t+5]==1)), values=c("dashed","solid","dotted")) + theme(legend.position=c(0.25, 0.75))
n4


#######################################################################################
### Figure S.12: Variations of long-term effects based on the counterfactual of interest for non-monotonic duration dependence (parabolic)
#######################################################################################
lte <- read.dta("LTE Example--NM1.dta", convert.underscore = TRUE)

### 1-time shift in effects
n1 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s1.0.0)) + xlab("") + ylab("Pr(Y=1)")
n1 <- n1 + geom_point() + geom_linerange(aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo), linetype="dashed")
n1 <- n1 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s1.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s1.1.0.hi, ymin = pr.s1.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n1 <- n1 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo, linetype="Y0"))
n1 <- n1 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n1 <- n1 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s1.1.0), vjust=-0.25, hjust=-.25, size=3)
n1 <- n1 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n1 <- n1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.75, 0.75))
n1

### Temporary (4-period) shift in effects
n2 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s2.0.0)) + xlab("") + ylab("Pr(Y=1)")
n2 <- n2 + geom_point() + geom_linerange(aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo), linetype="dashed")
n2 <- n2 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s2.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s2.1.0.hi, ymin = pr.s2.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n2 <- n2 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo, linetype="Y0"))
n2 <- n2 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n2 <- n2 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s2.1.0), vjust=-0.25, hjust=-.25, size=3)
n2 <- n2 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n2 <- n2 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.75, 0.75))
n2

### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.75, 0.75))
n3

### Compounded effects
n4 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n4 <- n4 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n4 <- n4 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n4 <- n4 + geom_point(data=subset(lte, xaxis > 5), aes(y=pr.s3.1.1, x=xaxis-0.3)) + geom_linerange(data=subset(lte, xaxis > 5), aes(ymax = pr.s3.1.1.hi, ymin = pr.s3.1.1.lo, x=xaxis-0.3, linetype="Y11"))
n4 <- n4 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.75, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 5), aes(label = time1.1, y=pr.s3.1.1), vjust=1, hjust=3.5, size=3)
n4 <- n4 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n4 <- n4 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1","Y11"), labels=c(expression(Y[t]==0), expression(Y[t]==1), expression(Y[t+5]==1)), values=c("dashed","solid","dotted")) + theme(legend.position=c(0.75, 0.75))
n4



#######################################################################################
### Figure S.13: Variations of long-term effects based on the counterfactual of interest for non-monotonic duration dependence (log-logistic)
#######################################################################################
lte <- read.dta("LTE Example--NM2.dta", convert.underscore = TRUE)

### 1-time shift in effects
n1 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s1.0.0)) + xlab("") + ylab("Pr(Y=1)")
n1 <- n1 + geom_point() + geom_linerange(aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo), linetype="dashed")
n1 <- n1 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s1.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s1.1.0.hi, ymin = pr.s1.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n1 <- n1 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s1.0.0.hi, ymin = pr.s1.0.0.lo, linetype="Y0"))
n1 <- n1 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n1 <- n1 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s1.1.0), vjust=-0.25, hjust=-.25, size=3)
n1 <- n1 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n1 <- n1 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.9, 0.9))
n1


### Temporary (4-period) shift in effects
n2 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s2.0.0)) + xlab("") + ylab("Pr(Y=1)")
n2 <- n2 + geom_point() + geom_linerange(aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo), linetype="dashed")
n2 <- n2 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s2.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s2.1.0.hi, ymin = pr.s2.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n2 <- n2 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s2.0.0.hi, ymin = pr.s2.0.0.lo, linetype="Y0"))
n2 <- n2 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n2 <- n2 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s2.1.0), vjust=-0.25, hjust=-.25, size=3)
n2 <- n2 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n2 <- n2 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.9, 0.9))
n2


### Permanent effects
n3 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n3 <- n3 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n3 <- n3 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n3 <- n3 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n3 <- n3 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n3 <- n3 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.25, hjust=-.25, size=3)
n3 <- n3 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n3 <- n3 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1"), labels=c(expression(Y[t]==0), expression(Y[t]==1)), values=c("dashed","solid")) + theme(legend.position=c(0.9, 0.9))
n3


### Compounded effects
n4 <- ggplot(subset(lte, xaxis >= 0 & xaxis <= 12), aes(x=xaxis, y=pr.s3.0.0)) + xlab("") + ylab("Pr(Y=1)")
n4 <- n4 + geom_point() + geom_linerange(aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo), linetype="dashed")
n4 <- n4 + geom_point(data=subset(lte, xaxis > 0), aes(y=pr.s3.1.0, x=xaxis-0.1)) + geom_linerange(data=subset(lte, xaxis > 0), aes(ymax = pr.s3.1.0.hi, ymin = pr.s3.1.0.lo, x=xaxis-0.1, linetype="Y1"))
n4 <- n4 + geom_point(data=subset(lte, xaxis > 5), aes(y=pr.s3.1.1, x=xaxis-0.3)) + geom_linerange(data=subset(lte, xaxis > 5), aes(ymax = pr.s3.1.1.hi, ymin = pr.s3.1.1.lo, x=xaxis-0.3, linetype="Y11"))
n4 <- n4 + geom_point(data=subset(lte, xaxis < 1)) + geom_linerange(data=subset(lte, xaxis < 1), aes(ymax = pr.s3.0.0.hi, ymin = pr.s3.0.0.lo, linetype="Y0"))
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(aes(label = time0.0), vjust=-0.25, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 0), aes(label = time1.0, y=pr.s3.1.0), vjust=-0.75, hjust=-0.25, size=3)
n4 <- n4 + geom_text(data=subset(lte, xaxis > 5), aes(label = time1.1, y=pr.s3.1.1), vjust=1, hjust=3.5, size=3)
n4 <- n4 + scale_x_continuous(limits=c(-2,12), breaks=c(-1, 0, 4, 8, 12), labels=c("Baseline", "t", "t+4", "t+8", "t+12"))
n4 <- n4 + scale_linetype_manual(name="Counterfactual", breaks=c("Y0","Y1","Y11"), labels=c(expression(Y[t]==0), expression(Y[t]==1), expression(Y[t+5]==1)), values=c("dashed","solid","dotted")) + theme(legend.position=c(0.1, 0.9))
n4



#######################################################################################
### Figure S.14: Average hazard rates compared to true hazard rates
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_1.dta", convert.underscore = TRUE)
mn_nonph <- subset(mn, scenario != "Negative with Positive NPH" & scenario != "Positive with Negative NPH")

p <- ggplot(subset(mn_nonph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85))
p

#######################################################################################
### Figure S.15: Average hazard rates compared to true hazard rates
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_2.dta", convert.underscore = TRUE)
mn_nonph <- subset(mn, scenario != "Negative with Positive NPH" & scenario != "Positive with Negative NPH")

p <- ggplot(subset(mn_nonph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85))
p


#######################################################################################
### Figure S.16: Average hazard rates compared to true hazard rates
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_3.dta", convert.underscore = TRUE)
mn_nonph <- subset(mn, scenario != "Negative with Positive NPH" & scenario != "Positive with Negative NPH")

p <- ggplot(subset(mn_nonph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85))
p


#######################################################################################
### Figure S.17: Average hazard rates compared to true hazard rates
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_4.dta", convert.underscore = TRUE)
mn_nonph <- subset(mn, scenario != "Negative with Positive NPH" & scenario != "Positive with Negative NPH")

p <- ggplot(subset(mn_nonph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85))
p


#######################################################################################
### Figure S.18: Estimated hazard rates compared to true hazard rates: negative duration dependence
#######################################################################################
hdata <- read.dta("Negative_Hazard_1.dta", convert.underscore=TRUE)
trueh <- subset(hdata, sim == 1)

### Temporal dummies
td <- ggplot(subset(hdata, t < 25), aes(x=t, y=TD.h, group=sim)) + ggtitle("Temporal Dummies") + ylab("Hazard") + xlab("")
td <- td + geom_line(color = "#d9d9d9")
td <- td + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
td

### Cubic polynomials
cp <- ggplot(subset(hdata, t < 25), aes(x=t, y=CP.h, group=sim)) + ggtitle("Cubic Polynomials") + ylab("Hazard") + xlab("")
cp <- cp + geom_line(color = "#d9d9d9")
cp <- cp + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
cp

### Cubic splines
s <- ggplot(subset(hdata, t < 25 & S.h != ""), aes(x=t, y=S.h, group=sim)) + ggtitle("Splines") + ylab("Hazard") + xlab("")
s <- s + geom_line(color = "#d9d9d9")
s <- s + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
s

### Automated smoothing splines
ass <- ggplot(subset(hdata, t < 25), aes(x=t, y=ASS.h, group=sim)) + ggtitle("Automated Splines") + ylab("Hazard") + xlab("")
ass <- ass + geom_line(color = "#d9d9d9")
ass <- ass + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
ass

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(td, vp = vplayout(1, 1))
print(cp, vp = vplayout(1, 2))
print(s, vp = vplayout(2, 1))
print(ass, vp = vplayout(2, 2))


#######################################################################################
### Figure S.19: Estimated hazard rates compared to true hazard rates: positive duration dependence
#######################################################################################
hdata <- read.dta("Positive_Hazard_1.dta", convert.underscore=TRUE)
trueh <- subset(hdata, sim == 1)

### Temporal dummies
td <- ggplot(subset(hdata, t < 25), aes(x=t, y=TD.h, group=sim)) + ggtitle("Temporal Dummies") + ylab("Hazard") + xlab("")
td <- td + geom_line(color = "#d9d9d9")
td <- td + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
td

### Cubic polynomials
cp <- ggplot(subset(hdata, t < 25), aes(x=t, y=CP.h, group=sim)) + ggtitle("Cubic Polynomials") + ylab("Hazard") + xlab("")
cp <- cp + geom_line(color = "#d9d9d9")
cp <- cp + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
cp

### Cubic splines
s <- ggplot(subset(hdata, t < 25 & S.h != ""), aes(x=t, y=S.h, group=sim)) + ggtitle("Splines") + ylab("Hazard") + xlab("")
s <- s + geom_line(color = "#d9d9d9")
s <- s + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
s

### Automated smoothing splines
ass <- ggplot(subset(hdata, t < 25), aes(x=t, y=ASS.h, group=sim)) + ggtitle("Automated Splines") + ylab("Hazard") + xlab("")
ass <- ass + geom_line(color = "#d9d9d9")
ass <- ass + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
ass

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(td, vp = vplayout(1, 1))
print(cp, vp = vplayout(1, 2))
print(s, vp = vplayout(2, 1))
print(ass, vp = vplayout(2, 2))


#######################################################################################
### Figure S.20: Estimated hazard rates compared to true hazard rates: non-monotonic duration dependence (parabolic)
#######################################################################################
hdata <- read.dta("NM1_Hazard_1.dta", convert.underscore=TRUE)
trueh <- subset(hdata, sim == 1)

### Temporal dummies
td <- ggplot(subset(hdata, t < 25), aes(x=t, y=TD.h, group=sim)) + ggtitle("Temporal Dummies") + ylab("Hazard") + xlab("")
td <- td + geom_line(color = "#d9d9d9")
td <- td + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
td

### Cubic polynomials
cp <- ggplot(subset(hdata, t < 25), aes(x=t, y=CP.h, group=sim)) + ggtitle("Cubic Polynomials") + ylab("Hazard") + xlab("")
cp <- cp + geom_line(color = "#d9d9d9")
cp <- cp + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
cp

### Cubic splines
s <- ggplot(subset(hdata, t < 25 & S.h != ""), aes(x=t, y=S.h, group=sim)) + ggtitle("Splines") + ylab("Hazard") + xlab("")
s <- s + geom_line(color = "#d9d9d9")
s <- s + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
s

### Automated smoothing splines
ass <- ggplot(subset(hdata, t < 25), aes(x=t, y=ASS.h, group=sim)) + ggtitle("Automated Splines") + ylab("Hazard") + xlab("")
ass <- ass + geom_line(color = "#d9d9d9")
ass <- ass + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
ass

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(td, vp = vplayout(1, 1))
print(cp, vp = vplayout(1, 2))
print(s, vp = vplayout(2, 1))
print(ass, vp = vplayout(2, 2))


#######################################################################################
### Figure S.21: Estimated hazard rates compared to true hazard rates: non-monotonic duration dependence (log-logistic)
#######################################################################################
hdata <- read.dta("NM2_Hazard_1.dta", convert.underscore=TRUE)
trueh <- subset(hdata, sim == 1)

### Temporal dummies
td <- ggplot(subset(hdata, t < 25), aes(x=t, y=TD.h, group=sim)) + ggtitle("Temporal Dummies") + ylab("Hazard") + xlab("")
td <- td + geom_line(color = "#d9d9d9")
td <- td + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
td

### Cubic polynomials
cp <- ggplot(subset(hdata, t < 25), aes(x=t, y=CP.h, group=sim)) + ggtitle("Cubic Polynomials") + ylab("Hazard") + xlab("")
cp <- cp + geom_line(color = "#d9d9d9")
cp <- cp + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
cp

### Cubic splines
s <- ggplot(subset(hdata, t < 25 & S.h != ""), aes(x=t, y=S.h, group=sim)) + ggtitle("Splines") + ylab("Hazard") + xlab("")
s <- s + geom_line(color = "#d9d9d9")
s <- s + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
s

### Automated smoothing splines
ass <- ggplot(subset(hdata, t < 25), aes(x=t, y=ASS.h, group=sim)) + ggtitle("Automated Splines") + ylab("Hazard") + xlab("")
ass <- ass + geom_line(color = "#d9d9d9")
ass <- ass + geom_line(data=subset(trueh, t < 25), aes(x=t, y=true.h), color="black")
ass

grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(td, vp = vplayout(1, 1))
print(cp, vp = vplayout(1, 2))
print(s, vp = vplayout(2, 1))
print(ass, vp = vplayout(2, 2))


#######################################################################################
### Figure S.22: Average hazard rates compared to true hazard rate under non-proportional hazards
#######################################################################################
mn <- read.dta("Average Hazard Rates--Scenario_3.dta", convert.underscore = TRUE)
mn_nph <- subset(mn, scenario == "Negative with Positive NPH" | scenario == "Positive with Negative NPH")

p <- ggplot(subset(mn_nph, t < 25), aes(x=t, y=mnh, group=funcform)) + xlab("Time") + ylab("Average Hazard")
p <- p + geom_line(aes(linetype=as.factor(funcform)), size=1)
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("1","2","3","4","5"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines","True Hazard"), values=c("dotdash","dotted","dashed","longdash","solid"))
p <- p + theme(legend.position=c(0.15,0.85))
p


#######################################################################################
### Figure S.23: Average predictive differences of X_K with alternative estimation techniques under conditions ignoring non-proportional hazards
#######################################################################################
nwp <- read.dta("NwP--APD Summary.dta", convert.underscore=TRUE)

n <- ggplot(subset(nwp, t < 20), aes(x=t)) + xlab("") + ylab("Average Predictive Differences")
n <- n + geom_line(aes(y=CP.APD, linetype="CP")) + geom_line(aes(y=TD.APD, linetype="TD")) + geom_line(aes(y=S.APD, linetype="S")) + geom_line(aes(y=ASS.APD, linetype="ASS"))
n <- n + facet_wrap( ~ scenario, ncol=2)
n <- n + scale_linetype_manual(name="Techniques", breaks=c("TD","CP","S","ASS"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines"), values=c("dotted","solid","dashed","longdash"))
n <- n + theme(legend.position=c(0.85,0.85))
n

#######################################################################################
### Figure S.24: Average predictive differences of X_K with alternative estimation techniques under conditions ignoring non-proportional hazards
#######################################################################################
pwn <- read.dta("PwN--APD Summary.dta", convert.underscore=TRUE)

p <- ggplot(pwn, aes(x=t)) + xlab("") + ylab("Average Predictive Differences")
p <- p + geom_line(aes(y=CP.APD, linetype="CP")) + geom_line(aes(y=TD.APD, linetype="TD")) + geom_line(aes(y=S.APD, linetype="S")) + geom_line(aes(y=ASS.APD, linetype="ASS"))
p <- p + facet_wrap( ~ scenario, ncol=2)
p <- p + scale_linetype_manual(name="Techniques", breaks=c("TD","CP","S","ASS"), labels=c("Time Dummies", "Cubic Polynomials", "Splines", "Smoothing Splines"), values=c("dotted","solid","dashed","longdash"))
p <- p + theme(legend.position=c(0.85,0.15))
p

#######################################################################################
### Figure S.25: Predicted probability of dispute initiation for coalition governments and ideological distance
#######################################################################################
rep = read.dta("Predicted Probabilities--Ideological Distance.dta", convert.underscore=TRUE) 
dist = read.dta("dist.dta", convert.underscore=TRUE) 

r <- ggplot(rep, aes(x = xaxis))
r <- r + geom_ribbon(aes(ymin = pr1.lo, ymax = pr1.hi), alpha = 0.3) + geom_line(aes(y = pr1))
r <- r + geom_line(aes(x = xaxis, y = pr0)) + geom_line(aes(x = xaxis, y = pr0.lo, type = 2), linetype = "dashed")  + geom_line(aes(x = xaxis, y = pr0.hi, type = 2), linetype = "dashed")
r <- r + geom_rug(data = dist, aes(x = dist.directional1), sides = "b")
r <- r + xlab("Ideological Distance: Outlier Party") + ylab("Pr(Initiation)") + theme(plot.title = element_text(family = "serif", size = 16, hjust = 0.5), axis.text.x = element_text(family = "serif", size = 12), axis.text.y = element_text(family = "serif", size = 12), axis.title.x = element_text(family = "serif", size = 14), axis.title.y = element_text(family = "serif", size = 14))
r

r2 <- ggplot(rep, aes(x = xaxis))
r2 <- r2 + geom_ribbon(aes(ymin = pr2.lo, ymax = pr2.hi), alpha = 0.3) + geom_line(aes(y = pr2))
r2 <- r2 + geom_line(aes(x = xaxis, y = pr0)) + geom_line(aes(x = xaxis, y = pr0.lo, type = 2), linetype = "dashed")  + geom_line(aes(x = xaxis, y = pr0.hi, type = 2), linetype = "dashed")
r2 <- r2 + geom_rug(data = dist, aes(x = dist.directional1), sides = "b")
r2 <- r2 + xlab("Ideological Distance: Outlier Party") + ylab("") + theme(plot.title = element_text(family = "serif", size = 16, hjust = 0.5), axis.text.x = element_text(family = "serif", size = 12), axis.text.y = element_text(family = "serif", size = 12), axis.title.x = element_text(family = "serif", size = 14), axis.title.y = element_text(family = "serif", size = 14))
r2

grid.newpage()
pushViewport(viewport(layout = grid.layout(1,2)))

vplayout <- function(x, y)
	viewport(layout.pos.row = x, layout.pos.col = y)
print(r, vp = vplayout(1, 1))
print(r2, vp = vplayout(1, 2))



#######################################################################################
### Figure S.26: Receiver-Operating Characteristic (ROC) Curve for Clare (2010)
#######################################################################################
sp = read.dta("Separation Plot.dta", convert.underscore=TRUE)

separationplot(sp$p, sp$init.monad, type = "line", line = TRUE)

#######################################################################################
### Figure S.27: Separation Plot for Clare (2010)
#######################################################################################
r = read.dta("ROC.dta", convert.underscore=TRUE)

roc <- ggplot() + geom_line(data = r, aes(x=prop1, y=prop0), colour = "blue") + geom_line(data = r, aes(x=x, y=y), linetype="dashed") + xlab("Proportion of 1s Correctly Predicted") + ylab("Proportion of 0s Correctly Predicted") 
roc


#######################################################################################
### Figure S.28: Box-whisker plots of the change in predicted probability for a 1-unit increase in X across the range of probabilities (0.2 range)
#######################################################################################
v = read.dta("Variability.dta", convert.underscore=TRUE)

v1 <- subset(v, s == 1)
v2 <- subset(v, s == 2)
v3 <- subset(v, s == 3)
v4 <- subset(v, s == 4)
v5 <- subset(v, s == 5)
v6 <- subset(v, s == 6)
v7 <- subset(v, s == 7)
v8 <- subset(v, s == 8)
v9 <- subset(v, s == 9)
v10 <- subset(v, s == 10)
v11 <- subset(v, s == 11)
v12 <- subset(v, s == 12)

g <- ggplot(v, aes(s, me)) 
g <- g + geom_boxplot(data = v1) + geom_boxplot(data = v2) + geom_boxplot(data = v3) + geom_boxplot(data = v4) + geom_boxplot(data = v5)
g <- g + ylab(expression(paste("Marginal effects of ", beta, "=1"))) + xlab("Probability Range") + coord_flip()
g <- g + scale_x_discrete(limits=c("[0.0, 0.2]","[0.2, 0.4]","[0.4, 0.6]","[0.6, 0.8]","[0.8, 1.0]"))
g


#######################################################################################
### Figure S.29: Box-whisker plots of the change in predicted probability for a 1-unit increase in X across the range of probabilities (0.4 range)
#######################################################################################
g2 <- ggplot(v, aes(s2, me)) 
g2 <- g2 + geom_boxplot(data = v6) + geom_boxplot(data = v7) + geom_boxplot(data = v8) + geom_boxplot(data = v9) + geom_boxplot(data = v10) + geom_boxplot(data = v11) + geom_boxplot(data = v12)
g2 <- g2 + ylab(expression(paste("Marginal effects of ", beta, "=1"))) + xlab("Probability Range") + coord_flip()
g2 <- g2 + scale_x_discrete(limits=c("[0.0, 0.4]","[0.1, 0.5]","[0.2, 0.6]","[0.3, 0.7]","[0.4, 0.8]", "[0.5, 0.9]", "[0.6, 1.0]"))
g2










